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The  extended  Capability  Magnus  Rotor  and  Ballistic  Body  6-DOF  Trajec¬ 
tory  Program  has  been  prepared  under  Contract  No.  F08635-70-C-0012  with  the 
Air  Force  Armament  Laboratory,  Eg! in  Air  Force  Base,  Florida,  by  Alpha 
Research,  Inc.,  Santa  Barbara,  California.  The  programmer  at  Alpha 
Research,  Inc.  was  Mr.  William  Davidson.  The  program  monitor  for  the 
Armament  Laboratory  was  Mr.  Edward  S.  Sears  (ATRA).  This  effort  was  con¬ 
ducted  during  the  period  6  October  1969  to  6  April  1970. 

This  program  is  a  modification  of  several  earlier  computer  programs 
prepared  for  magnus-rotor  flight  dynamics  investigations.  The  original 
program  was  prepared  for  the  U.  S.  Army  Edgewood  Arsenal  under  Contract  No. 
DA-1 8-1 08-AMC-236(A) .  The  original  program  was  later  adapted  to  the  Air 
Force  Armament  Laboratory  computer  facilities  by  Air  Force  personnel,  and 
was  used  by  Alpha  Research,  Inc.  in  conjunction  with  Air  Force  Contracts 
Nos.  F08635-67-C-01 35  and  F08635-69-C-0106. 

The  present  modified  program  is  written  in  General  Fortran  IV  language 
compatible  with  third  generation  computers  such  as  the  GE  635,  IBM  360, 
and  CDC  6400. 

Information  in  this  report  is  embargoed  under  the  Department  of  State 
International  Traffic  In  Arms  Regulations.  This  report  may  be  released  to 
foreign  governments  by  departments  or  agencies  of  the  U.  S.  Government 
subject  to  approval  of  the  Air  Force  Armament  Laboratory  (ATRA),  Eglin  AFB, 
Florida  32542,  or  higher  authority  within  the  Department  of  the  Air  Force. 
Private  individuals  or  firms  require  a  Department  of  State  export  license. 

This  technical  report  has  been  reviewed  and  is  approved. 

CHARLES  K.  ARPKE,  Lt  Colonel,  USAF 
Chief,  Technology  Division 
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ABSTRACT 


A  six-degrees-of-freedom  trajectory  program  for  quasi -symmetric 
rigid  bodies  is  described.  The  equations  of  motion  are  developed  such  that 
either  a  body-fixed  or  fixed-plane  moving  coordinate  system  can  be  utilized. 
Provision  is  made  for  large  angle  and  angular  rate  motions,  such  as  are 
experienced  by  magnus  rotor  munitions. 

The  aerodynamic  system  permits  the  usual  aeroballistic  coefficients 
to  be  expressed  as  tabular  functions  of  angle  of  attack  and  Mach  number; 
in  addition^the  magnus  force,  magnus  moment,  and  rolling  moment  coeffi¬ 
cients  can  be  tabular  functions  of  the  nondimensional  spin  parameter,  pd/2V. 
Additional  aerodynamic  terms  are  provided  to  account  for  body-fixed  aero¬ 
dynamic  asymmetries  and/or  control  inputs,  aerodynamic  roll  angle  effects, 
flow  asymmetry  with  respect  to  the  angle  of  attack  plane  at  zero  spin,  and 
lateral  c.  g.  offset. 

The  computer  program  is  written  in  General  Fortran  IV  language 
compatible  with  CDC  6400,  IBM  360,  and  GE  635  data  processing  machines. 
Included  in  the  report  are  the  program  input  and  output  formats,  flow  charts 
of  the  main  program  subroutines,  and  a  complete  program  listing. 
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SECTION  I 


PROGRAM  DESCRIPTION 


A.  INTRODUCTION 

The  present  computer  program  has  resulted  from  the  need 
for  more  exact  simulation  of  the  motion  of  dispenser  munitions.  Particular 
attention  has  been  given  the  simulation  of  the  flight  of  magnus  rotor  bomb- 
lets.  The  program  can  also  provide  trajectory  and  motion  simulation  for 
most  unpowered  projectiles,  missiles,  and  rockets. 

The  most  significant  features  of  the  extended  capability 
trajectory  program  are: 

1.  Choice  of  fixed-plane  or  body-fixed  axes. 

2.  All-attitude  motion  prediction. 

3.  Adaptability  to  very  large  spin  rates. 

4.  All  basic  aeroballistic  coefficients  are  tabular 
functions  of  o2  and  M .  In  addition,  spin  and  magnus 
coefficients  are  tabular  functions  of  pd/2V. 

5.  Inclusion  of  high  order  nonlinear  damping  terms. 

6.  Aerodynamic  dependence  upon  roll  angle  is  included. 

7.  Provision  for  aerodynamic  and  configurational 
asymmetries,  including  c.  g.  lateral  offset. 

8.  Provision  for  wind. 

9.  Provision  for  elimination  of  high  frequency 
(nutational)  motion. 


B.  COORDINATE  SYSTEMS 

Either  of  two  moving  axes  systems  can  be  selected,  a  body- 
fixed  coordinate  system  or  a  fixed-plane  coordinate  system.  Both  sets  of 
moving  coordinates  have  as  their  origin  the  center  of  mass  of  the  body. 
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The  inertial  reference  system  is  a  set  of  XYZ  non-rotating 
earth-fixed  coordinates,  with  the  origin  at  sea  level.  The  orientation  of 
the  moving  axes  with  respect  to  the  inertial  axes  is  given  by  three  Euler 
angles:  Q  ,  '\p~  ,  and  (f>  ,  as  depicted  in  Figure  1, 

Body-Fixed  Coordinates  The  body -fixed  coordinates  are 
comprised  of  the  right-handed  orthogonal  set  of  axes,  x  ,  y  ,  and  z  . 

Fixed-Plane  Coordinates  The  fixed-plane  coordinates  are 
comprised  of  the  right-handed  orthogonal  set  of  axes  x  ,  y'  ,  and  z'  .  In  the 
fixed-plane  system,  the  y1  axis  is  initially  in  the  XY  plane  and  is  so  re¬ 
strained  as  to  stay  in  that  plane.  Consequently,  (p  =  q>  »  O  for  the  fixed- 
plane  coordinate  system. 

Since  the  body  can  rotate  with  respect  to  the  x  y'  z'  axes, 
this  axes  system  is  restricted  to  bodies  with  rotational  symmetry. 


C.  ATTITUDE  REPRESENTATION 

The  orientation  of  the  moving  coordinates,  for  input  and 
output  purposes,  is  described  by  the  Euler  angles.  For  computational  pur¬ 
poses,  however,  an  angular  orientation  scheme  based  on  quaternions  is 
used.  The  quaternion  system  avoids  the  discontinuities  which  occur  in  the 
trigonometric  functions  and  angular  rates  when  the  motion  passes  from 
quadrant  to  quadrant. 

The  scalar  relations  between  the  four  parameter  quaternion 
system  and  the  three  Euler  angles  are  defined  for  the  body -fixed  and  fixed- 
plane  system  as  follows:^) 

Body -fixed  axes: 

%)  C^(e/ *)««.(*&)  +  ** 

A,  -  AW  (#L  (*%.)  CM 0%)  -  CM  (%)  Awfflz)  MmWz) 
X2  -  (*/* )  AM  (*/2)  4-  AM  (%)  U»k)  Aw  (Vt) 

/l  3  “  CM  (f/z)  C&4  W  AW  C%)  ~  Aw  m  Aw  Wt)  CM  IVz) 
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:  x  is  body  longitudinal  axis 

or  axis  of  rotational 
symmetry 


z 


Inertial  Axes 


Figure  1.  Coordinate  Systems 


Fixed- plane  axes: 


A0  =  ocr^  (Q/z)  cc*. (Wz) 

A,  =  -  Mm*  {g/z  )  Ju*  [T/z  ) 

A2  *  AOk  (&/2)  Ctr^Wz) 

X3  =  {®!z) Mm.(VVz) 

Euler  angles : 

Mai  0  =  2  f  A  A2-A,A3)  ’  -  ft/z  -  0  ^  K/z 

hbK,  6  - 

U*Y  =  t A^As)  -  ~n  <  v  <>  ir 

A/+Af  -A, -A3 

to* 1  <T)  =  _2.CAa.A3  4-  Ao  A,  )  ,  „<j-  ^  Q)  <  fl- 

A*  "  A,1  -A*  +  A|  '  * 

Rotation  Matrices  In  the  equations  of  motion,  transforma¬ 
tions  from  the  moving  axes  to  the  fixed  inertial  axes  are  required.  These 
can  be  expressed  (for  either  axis  system)  by  the  general  quaternion  rotation 
matrix  given  below,  where  the  subscripts  F  and  M  refer  to  the  fixed 
inertial  axes  and  the  moving  axis,  respectively. 


0 

/  0  00 

0 

?F, 

0  A.+\— A.-A,  2fx(A, -A. A,)  2  CA,  A,  ■*'A»A,1 

5f, 

O  2(A,Aa+\.A,A  X^^-A,1  2(A,A,-A.A,) 

.  V 

0  2(A,A3-A„A2)  2 (a, A,  *■  A, A,)  X.’-A’-A^ 

1 

C/rt 

Z 

_ 

_ 2  ( Ao  A2  ~~  Ai  A3) _ 
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D. 


BASIC  EQUATIONS  OF  MOTION 


The  equations  of  motion  consist  of  13  differential  equations 
for  the  variables 


Equations  of  Motion 


which  must  be  integrated  to  obtain  the  desired  motion  solution  and  trajectory. 
The  variables  u,  v,  and  w  are  the  linear  velocity  components  in  the  direc¬ 
tion  of  the  moving  axes,  and  p  ,  q  ,  and  r  are  the  angular  velocity  com- 
ponenty  corresponding  to  the  moving  axes.  Following  standard  notation, 
these  components  will  be  with  respect  to  the  x  ,  y  ,  z  or  x  ,  y1  ,  z1,  depend¬ 
ing  upon  whether  the  moving  axes  are  body -fixed  or  fixed-plane,  respectively. 
It  follows  that  p  and  u  are  identical  in  the  two  coordinate  systems,  but 
v  ,  w  ,  q  ,  and  r  are  not. 

The  derivation  of  the  fixed-plane  equations  of  motion  must, 
of  necessity,  consider  the  simplified  inertial  tensor 

■  Mf.p. 


Ix  o  0 

0  I  0 
0  0  I 
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The  body -fixed  equations  of  motion  are  more  general. 
Provision  is  made  for  bodies  with  an  inertia  tensor  * 


■'•x  *xy  0 


Wy  0  “  Hb.it. 


0  0 


corresponding  to  the  xy  plane  as  a  plane  of  mirror  symmetry.  This  is 
sufficient  for  most  simulation  work. 

From  basic  mechanics,  we  let  XI  represent  the  angular 
velocity  vector  of  the  coordinate  system  with  respect  to  the  inertial  system 
and  CO  represent  the  angular  velocity  of  the  body  with  respect  to  the  moving 
coordinates,  and  the  fundamental  equations  of  motion  become  (for  constant 
mass  and  inertia) 


where 


F  =  -mg  +  mV  +  S\  X  mV 
M  =  I  to  +  il  X  [l]<*> 


ZJ  B.F. 


a'J  F.P. 


=  aerodynamic  force 


M  = 


B.F. 


F.  P. 


=  aerodynamic  moment 


j  B.  F. 


fz’J  F.P. 


=  gravitational  acceleration 


B.  F. 


F.  P. 


=  velocity 


All  products  of  inertia  are  positive  quantities,  I„  =  I  , 

*y  yx 
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GO 


•  — 

p 

r  tan  Q 

q 

= 

q 

_  r 

B.  F. 

L  r  . 

P 

P 

q 

= 

q 

r 

B.  F. 

r 

angular  velocity  of 
coordinates 


angular  velocity  of  body 


and  tan  Q 


K-K-K+K 


Substitution  of  the  appropriate  quantities  into  the  fundamental  equations  of 
motion  results  in  the  scalar  equations  for  u,  v.  w,  p,  q,  and  r  .  The 
differential  equations  for  X  ,  Y  ,  Z  are  obtained  by  transformation  of  the 
components  u  ,  v  ,  w  .  Finally,  the  derivatives  of  the  quaternions  are 
obtained  from  ( 

A  -  i  A  *il 


where  *  denotes  a  non-commutative  quaternion  product, 
a*b  =  (aQ  +  aji  +  a2j  +  a3k)  (bQ  +  bji  +  b2j  +  b3k) 
which  can  be  expressed  in  matrix  form  as 


<a0)  (-a^)  (-a2)  (-a3) 

’bo" 

(a^)  (aQ)  (-a3)  (  a2 ) 

bl 

(a2)  (a3)  (aQ)  (-a^ 

-  b2 

(a3)  (-a2)  {  aJ  )  (  aQ  ) 

_  b3_ 

The  resulting  equations  of  motion  are  summarized  in  Tables 
I  and  II  for  the  body -fixed  and  fixed-plane  axes,  respectively. 

Note  that  a  singularity  exists  in  the  fixed-plane  equations  of 
motion  for  9  =  t  It /2,  which  precludes  the  use  of  these  equations  of 
motion  when  9  may  approach  7T/2. 
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TABLE  I.  BODY- FIXED  EQUATIONS  OF  MOTION 


• 

tX 

=  f  IT 

-  fur  +  2  f A,A3 

-  A*  Ax )  g  *  Fx  />vx 

r 

'  f“" 

-  ru  +  2  f  A2  A3 

+  A*A,)9  +  Fy/* 

Ar 

+  (X-  A,v 

"  A*  +  A3)g  +  Fj/m 

^  (I* 

r)-(^]r  + 

L 

lx 

.  -Kr« 

+  ^  “  Jsl  3xy  1  P  +  [ 

Ixy  +  4’?*  "*•  -ly'L  -*•  In^  H 

lxl_y  ~  Xcy* 

• 

r 

^fr 

+*rHvy 

+  £- 
TJ 

r(r»+i>-r3V»rr]fr-  [  V  *  rx(rx-i,)>-r  4VL  -h,«h 

lx  Iy  -  -I xy 

r  =  +  M 

X»  (x.l+  xr  -  xt  -  As  )  a  +  2(A,Az-AeA3)tr  +  z(rA1A3 +A0A,W 

V  *  2(a,A1  4>0As)  lc  +  (A*-A* +£- A,)ir  4- 2  (A*  A., -A.  A, V* 
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TABLE  II.  FIXED- PLANE  EQUATIONS  OF  MOTION 


CL 

ir 

Ur 


fw  -  for  +•  2  (  A,  As  -UOs  +  R  A 

2  UT  "  XAi  _  f.t4.  +  In. 

r-^r  w\i-  ^  ^  3 + 4 


•  — 
7  I, 


7 


„  f  A, A3  -  A. Ax  _  Ix^  M 

■  T  ir'  T 


-  f  (  2-1C 


A,  A^  —  X  Aj. 

A)- A,1-  Ax  +A3 


I  '  /  I 

rV4 


1* 

I 


X  = 

y  - 
z  - 

X.- 
A,  0 
V- 


( A.1  +  X- Ai-Aj)*  +  2  (A(At- A.A*)^  +  2 ( A.Aj  +  >.A1)ur 
zfA.Ax  +  XX)  «  4  (  A.'- AC"  +  A^-AsV  4  2(A,A3- A.A,)ar 
2  (  A, A3  -  A*)0«  4  2(  AaA3  +  A.  A,V  4  (A*-  A'- A^+AjV 

"  ^  (  M  *  X-X-X 
~  ^  *  a:-a,m:+a^  ) 

^  +  >1- A" -  x  +x  ) 


9 


E. 


BASIC  AERODYNAMIC  SYSTEM 


Basic  Aerob&ilistic  Coefficients  The  aerodynamic  forces 
and  moments  are  expressed  in  coefficient  form,  using  an  aeroballistic  sysw* 
tern  consistent  with  that  which  is  used  for  symmetric  missiles,  (^  The 
basic  aeroballistic  coefficients  are  valid  for  either  the  body -fixed  or  fixed- 
plane  coordinate  systems,  and  only  require  that  the  configuration  have 
trigonal  or  greater  rotational  symmetry.  These  basic  aeroballistic  coeffi- 


axial  force  coefficient 
normal  force  coefficient 
overturning  moment  coefficient 


cients  are  given  below: 

Cx 

(  Si  , 

M) 

CN 

(  Si  , 

M) 

CM 

(&  • 

M) 

CMq 

(Si  , 

M) 

CNp 

(  Si  , 

M, 

£*) 
2V ' 

CMp 

(  &L  . 

M, 

2V* 

CJ ! 

(«•  , 

M, 

2Y  ' 

C*p 

(  Sc  , 

M, 

El) 
2V ' 

=  pitch  damping  coefficient  based  on 

=  “^Lp  (  ex.  ,  M,  2^)  =  magnus  force  coefficient 
for  body -fixed  and  fixed-plane  axes,  respectively. 

=  magnus  moment  coefficient 

=  spin  torque  coefficient  due  to  canted  fins  or  vanes 


The  above  coefficients  depend  only  upon  the  total  angle  of 
attack,  oC  ,  and  are  independent  of  the  components  of  the  angle  of  attack.  * 
These  coefficients  are  also  functions  of  Mach  number,  and  in  addition,  Cjj  , 
CMp,  Cj ,  and  Cjp  are  permitted  to  be  functions  of  the  nondimensiona!  ” 

spin  parameters  ■  pd/2V. 

The  sense  of  the  basic  aeroballistic  forces  and  moments  is 
indicated  in  Figure  2.  The  aeroballistic  forces  and  moments  are  trans¬ 
formed  to  the  forces  and  moments  corresponding  to  the  moving  coordinates 


*  The  total  angle  of  attack  is  the  angle  between  x  axis  and  the  aerodynamic 
velocity  vector. 
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by  the  matrices 


"Fx 

V 

"  1  0  0 

’cx 

Fy 

sr 

Fy' 

= 

0  {-cos  S  )  (si*1  £  ) 

CN 

F*. 

.f*L 

0  {-sin  |  )  (-cos  %  ) 

CN 

L  NP-I 

j  p  V“S 


L 

•  — 

1  0  0 

Cl 

M 

x 

0  (sin  |  )  (cos  S  ) 

CM 

N 

0  {-cos  g )  (sin  S  ) 

.  CMp  . 

■  f  p  y2sd 


where 


cos 


I  = 


AT 


4at2+ 


our * 


sin  |  =  - 


ur 


Jatz+  ur* 


Additional  Damping  Parameters  In  general,  the  cross - 
angular  velocity  does  not  coincide  with  the  plane  of  the  total  angle  of  attack. 
For  these  nonplanar  motions,  it  is  now  generally  accepted  that  the  aero¬ 
dynamic  damping  can  differ  from  that  for  a  planar  motion.  ^  The  present 
trajectory  program  accounts  for  the  nonplanar  motion  damping  by  dividing 
the  cross -angular  velocity  into  components  q'  and  r',  which  are  coincident 
with,  and  normal  tg0the  angle  of  attack  plane,  respectively.  The  nonplanar 
damping  contribution  due  to  r1  is  incorporated  by  use  of  an  additional  damp¬ 
ing  coefficient,  Cnr.  In  a  similar  manner,  the  coupling  coefficients,  Cm, 


and  Cnp(j» 
Figure  3. 


are  introduced.  The  sense  of  these  coefficients  is  depicted  in 


pr 


The  moments  corresponding  to  the  moving  axes  are  deter¬ 
mined  by  the  double  rotation  matrices? 


M 

(sin  %  )  (cos  S ) 

(sin  g  )  (-cos  |) 

N 

{-cos  %  )  (sin 

(cos  g)  (sin  g  ) 

.<#>•  Cnr  . 

M 

(sin  g  )  (cos  %  ) 

(sin  g)  (cos  g  ) 

(fl*  '  CmPr 

N 

(-cos  g  )  (sing) 

(-cos  g )  (sin  g) 

1  2 
2  f°V Sd 

£  py2sd 
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are 


The  damping  coefficients,  Cnr ,  Cmpr  >  and  Cnp^  , 
functions  of  total  angle  of  attack  and  Mach  number,  as  summarized  below: 


Cnr  (  c*.  ,  M)  =  non  planar  damping  coefficient 


Cmpr  (  oc  ,  M)  =  cross  damping  moment  -  angle  of  attack  plane 


Cnpq  (£.  M) 


cross  damping  moment  -  magnus  plane. 


These  coefficients,  together  with  the  basic  aeroballistic 
coefficients,  comprise  the  aerodynamic  coefficient  system  for  the  fixed- 
plane -axes  system. 


F.  EXTENDED  AERODYNAMIC  SYSTEM 

When  the  body -fixed -axes  option  is  selected,  additional 
aerodynamic  coefficients  are  included  to  account  for  body-fixed  asymme¬ 
tries,  windward  meridian  orientation  (aerodynamic  roll  angle),  flow 
asymmetry,  and  lateral  displacement  of  the  center  of  gravity.  These 
coefficients  are  depicted  schematically  in  Figure  4. 

Body-Fixed  Asymmetries  Body-fixed  aerodynamic  forces 
and  moments  due  to  misalignment,  cant,  or  asymmetry  of  body  and/or 
lifting  surfaces  are  accounted  for  by  t'  e  coefficients 

Cyo  (£.  M) 

CZo  (Si,  M) 

^m0  ( ot  »  M) 

CnG  (  «*■ »  M>  * 

The  s.erodyn?.mic  effects  of  misalignment,  cant,  and  asymme¬ 
try  on  the  axial  force  and  rolling  moment  are  accounted  for  by  the  basic 
aeroballistic  coefficients  Cx  (3.  ,  M)  and  (  ex. ,  M). 

Aerodynamic  Effects  Due  to  Windward  Meridian  Orientation. 
Bodies  with  wings  and  fins  are  subject  to  aerodynamic  effects  related  to  the 
orientation  of  the  aerodynamic  surfaces  with  respect  *o  the  windward  merid¬ 
ian  of  the  cross  flow.  The  orientation  of  the  aerodynamic  surfaces  with 
respect  to  the  cross  flow  is  specified  by  the  aerodynamic  roll  angle  $  . 

This  angle  is  defined  as  a  clockwise  rotation  of  the  aerodynamic  surfaces 
with  respect  to  the  cross  flow,  looking  in  the  direction  of  the  x  axis,  as 
described  by  the  following  sketch.  The  orientation  of  the  symmetry  planes 
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Figure  4.  Additional  Aerodynamic  Forces  and  Moments 
for  the  Body- Fixed  Axes  Option 


of  the  fins  or  wings  with  respect  to  arbitrary  body  axes  is  given  by  the  phase 
angle  J  .  This  phase  angle  may  be  necessary  if  body-fixed  asymmetries, 
such  as  lateral  c.  g.  offset,  are  also  present. 


The  most  significant  induced  aerodynamic  effects  resulting 
from  windward  meridian  are  generally  the  induced  rolling  moment,  the 
induced  side  force,  and  the  induced  side  moment.  The  two  moments  often 
have  a  profound  effect  on  vehicle  stability,  and  are  usually  of  greatest  mag¬ 
nitude  for  3  or  4- fin  tail  assemblies.  Likewise,  very  large  induced 
aerodynamic  effects  would  be  expected  from  a  plan:  r  surface,  such  as  a 
wing.  In  addition,  the  aerodynamic  roll  angle  will  have  some  effect  on  the 
normal  force  and  pitching  moment,  which  must  be  considered  for  complete¬ 
ness. 

For  a  single  set  of  aerodynamic  surfaces  of  uniform  size,  the 
functional  dependence  of  the  aerodynamic  forces  and  moments  on  the  aero¬ 
dynamic  roll  angle  is  given  by  the  following  coefficients  and  functional 
relationships : 


lb 


CSFi  (  c*  .  M,  $  )  =  Csfx  (  <*.  ,  M)  sin  (S'),  , 

(  «£  ,  M,  $  )  =  Cjvj1  (  «£  ,  M)  sin  ^ 

CSMi  (  <*  »  M>  $  )  =  CSMj  (  S  ,  M)  sin 

CMi  (  OL  »  M»  $  )  =  (  CL  ,  M)  sin  (t^  ,  $  ^ 

Cjj(  (  SL  ,  M,  $  )  =  Ci|(  (  5.  ,  M)  sin  (V\1  ,  $  ^ 

where  V\  =  number  of  axially  symmetric  fins  (i.  e.  ,  cruciform  fins,  \  =  4). 

Although  the  dependence  of  CgFl>  CNl>  CgMl>  CMl»  and 
C#_  on  $  could  be  more  general,  only  the  first  harmonic  is  used.  The 

principal  reasons  for  this  simplification  is  the  usual  lack  of  more  definitive 
wind  tunnel  data  (i.  e. ,  usual  practice  is  to  measure  the  near  maximum 
effect  of  $  by  testing  at  $  =  IT  /2J[). 

Limited  provision  is  made  for  a  second  set  of  aerodynamic 
surfaces  with  a  different  >1  .  For  the  second  set  of  aerodynamic  surfaces 
only  the  coefficient  is  provided. 

*•2 

Experience  has  shown  that,  in  general,  one  set  of  aerodynamic 
surfaces  will  have  aerodynamic  induced  effects  which  are  much  larger,  and 
the  complete  force  and  moment  coefficients  should  be  used  with  this  set  of 
aerodynamic  surfaces. 


Combined  Effect  of  Geometric  Asymmetry  and  Windward 
Meridian  Orientation  The  preceding  paragraphs  describe  how  provision 
has  been  made  for  both  the  body-fixed  aerodynamic  coefficients  CyQ,  C Zq , 
Cm0i  and  Cno  and  windward  meridian  orientation.  In  practice,  a  single 
vehicle  modification  (or  asymmetry)  may  lead  to  both  effects  simultaneously. 
And  it  is  important  to  recognize,  clearly,  how  the  aerodynamic  coefficients 
should  be  interpreted. 


Consider  a  body  with  a  portion  of  the  nose  sliced  off  as  in 
Figure  5,  such  that  the  xz  plane  is  a  plane  of  mirror  symmetry.  Then,  in 
addition  to  the  body- fixed  trim  force  CZq  (  aZ.  )  and  the  normal  force 
Cn  (  o£  )  it  is  possible  that  induced  aerodynamic  normal  and  side  forces 
exist.  These  are  Cj^  (  <5.  ,  5  )  and  CgF^  (  oc  ,  $  ),  respectively,  and 

have  the  vector  orientations  indicated  by  Figure  5.  Note  further  that  in 
Figure  5,  which  depicts  a  case  for  which  \  =  1,  the  first  harmonic 
relationship  used  for  Cj<ji  and  CgFj  requires  that  Cj^j  =  C.F.  =  0  at 

§  =  9>,  I**  •.  , 
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Aerodynamic  Cross  Force  Due  to  Combined  Effect  of 
Geometric  Asymmetry  and  Windward  Meridian  Orientation 


if-'p 


At  this  point,  it  is  very  important  to  observe  that  the  effect  of 
the  body- fixed  trim  force  may  be  incorrectly  interpreted  as  having  a  $ 
dependence.  For  example,  if  the  y  z  axes  are  fixed  and  there  is  a  rotation 
of  the  windward  meridian,  then  there  ill  occur  a  corresponding  angular 
rotation  between  the  vectors  C^j  (  Si  )  and  CZq  (  Si  ),  such  that  the  magni¬ 
tude  of  the  total  cross  force  will  vary  with  §  .  This  effect  is  only  due  to 
the  variation  of  the  normal  force  with  §  ,  and  is  correctly  accounted  for  in 
the  body- fixed  equations  of  motion,  due  to  the  fact  that  FycC  -Cjvj(  Si  )  cos  £ 
and  Fz  oC  =  Cjvj  (  oL  )  sin  J  .  Thus,  in  reducing  wind  tunnel  force  measure- 
ments  to  determine  the  effect  of  $  ,  when  configurational  asymmetries  are 
present,  the  trim  force  coefficients  CyQ  and  CZq  must  be  subtracted  first. 
Similar  arguments  apply  to  the  aerodynamic  moment  coefficients. 

Additional  Aeroballistic  Coefficients  for  Flow  Asymmetry 
The  aerodynamic  forces  on  a  pure  axi- symmetric  body  (i.  e. ,  body  of  revo¬ 
lution)  with  zero  spin,  are  not  always  symmetric  with  respect  to  the  angle 
of  attack  plane  as  postulated  in  the  basic  aeroballistic  formulation.  In 
particular,  it  has  been  noted  that  long  pointed-nose  bodies  tend  to  have  an 
asymmetric  vortex-wake  structure  at  large  angles  of  attack  for  a  range  of 
Reynolds  numbers  and  Mach  numbers.  Once  established,  the  asymmetrical 
wake  becomes  reasonably  stable,  such  that  the  resulting  forces  and  moments 
can  be  considered  steady.  To  include  the  above  effects  in  the  aerodynamic 
model,  the  following  coefficients  and  functional  relationships  are  added  to 
the  aeroballistic  system  employed  with  the  body- fixed  axes: 

CSF3  (  oc  ,  M  ) 

CSM3  (  5.  f  M  ) 

It  will  be  noted  that  these  coefficients  have  the  same  vector 
orientation  as  the  magnus  force  and  moment. 

Effect  of  Lateral  Displacement  of  Center  of  Gravity  from  the 
Longitudinal  Reference  Axis  Additional  aerodynamic  moment  terms  are 
introduced  into  the  equations  of  motion  if  the  center  of  gravity  is  laterally 
displaced,  from  the  longitudinal  reference  axis. 

These  moments  are: 

AC*  =  -(  E  Cz)  (  Ay/d) 

ACn  =  Cx  (  Sc.  ,  M)  (  Ay/d) 

Only  a  c.  g.  offset  along  the  y  body- fixed  axis  will  be  con¬ 
sidered,  because  all  other  body-fixed  forces  and  moments  are  introduced 
with  complete  generality.  Also,  the  c.g.  offset  conforms  to  the  xy  plane, 
which  is  assumed  to  be  a  plane  of  mirror  symmetry. 
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G. 


SUMMARY  OF  AERODYNAMIC  FORCES  AND  MOMENTS 


Table  III  summarizes  the  aerodynamic  forces  and  moments, 
in  scalar  form,  for  each  degree  of  freedom.  Those  terms  which  are  added 
for  the  body-fixed  axes  option  are  enclosed  by  dashed  lines. 


H.  ATMOSPHERIC  MODEL  AND  WIND  SIMULATION 


Atmospheric  Model  Air  density,  p  ,  and  velocity  of  sound, 
a  ,  are  approximated  for  standard  day  conditions  by  the  relations: 


P 

P 

a 

a 


0.0023769  [  1  +  6.875  x  1 0-6 (Z) ]4' 2561 


0.0040  e 


4, 806  x  10"D(Z) 


968.46  -  0.  004123  (:Z) 
968.46 


(-Z)  <  36,  000  ft. 

(-Z)  >  36,  000  ft. 

(-Z)  <  36,  000  ft. 
(-Z)  >  36,  000  ft. 


I.  ATMOSPHERIC  WIND  SIMULATION 

Atmospheric  wind  simulation  is  accomplished  by  introducing 
wind  vectors  Xw  ,  Yw  ,  as  a  function  of  altitude  (  -Z).  The  aerodynamic 
velocity  components  of  the  body  with  respect  to  the  moving  air  mass  are 


(X  -  Xw)  ,  (Y  -  Yw)  ,  Z 


and  the  corresponding  aerodynamic  velocity  components  along  the  body-axes 
are  designated  u^  ,  v^  ,  w^  .  In  the  presence  of  wind,  the  aerodynamic 
forces  and  moments  are  redefined  so  as  to  be  functions  of 


2 


M 


-*  Mh 

A 

a 


I 


tiM.  1 


UJa 
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where 


Mm.  a  ~ 


Ctr*-  t  = 


Ak 


The  above  relations  are  of  the  same  algebraic  form  as  in  the 
case  of  zero  wind. 

The  transformation  from  the  inertial  wind  aerodynamic 
velocities  to  the  body-axes  aerodynamic  velocities  is  given  by 


2(A,Ai  +  AbA3) 

2.(A,  X3  -A*  A,) 

= 

2.(  A,Aj  ~  A,A3} 

A*  -  A,  +  A2  -  A3 

%•(  A*  A3  +  A0A,) 

"a 

•  «* 

^(A,A3  Ab  A*) 

2.  C  A2  A3  -  A.  A,) 

a:-a:-a:+x23 

X  *  x« 
Y-Vw 
Z*  Zw 


J.  QUATERNION  NORMALIZATION  AND  QUATERNION  ERROR 

The  separate  integration  of  each  quaternion  element,  X,“^Aj, 
leads  to  the  possibility  of  small  errors,  such  that  X„  +  A*  +  +  X3  I  . 

To  insure  normalization  (  X*  +  A*  +  X*  +•  Xj  =  I  )  ,  the  following 
correction  is  made  to  each  quaternion  element  after  integration: 


The  above  correction  scheme  is  derived  from  the  transforma¬ 
tion  relationships  between  the  quaternions  and  the  three  Euler  angles,  which 
uniquely  describe  the  moving-axis-system  orientation.  The  method  results 
in  an  exact  normalized  quaternion.  However,  in  the  case  of  the  fixed-plane 
axes,  tue  additional  identity,  X0X,  +  Xj  A3  =S  0  ,  which  applies  for  this  case, 
may  not  be  satisfied.  Consequently,  the  body-fixed  axes  should  be  used  in 
preference  to  the  fixed-plane  axes  whenever  the  rotational  rates  and  frequen¬ 
cies  permit. 
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The  quaternion  error.  £  ,  is  defined  as  the  vector  sum  of 

the  errors  of  the  normalized  quaternion  elements  after  successive  predict- 
correct  or  correct  and  re-correct  integration  operations:  i.  e. , 


£ 


2 


normalized  value  of  A^  after  prediction  (or  correction; 
normalized  value  of  A.j  after  correction  (or  re- correction) 

K.  INTEGRATION 

Since  the  differential  equations  are  quite  complicated  and 
lengthy,  a  refined  integration  scheme  is  employed.  The  basic  system 
utilizes  Milne^  four-point  method  of  prediction  and  Simpson's  rule  for 
correction.  The  Runge-Kutta  method  of  third  order  accuracy  is  used  to 
calculate  the  second  through  fourth  ordinates  needed  to  start  Milne's  method. 
These  are  described  in  detail  in  the  flow  charts.  Section  IV. 

Options  are  also  provided  for  the  Adams  and  trapezoidal 
integration  schemes,  using  four  and  three  ordinates,  respectively. 

As  a  further  means  of  reducing  integration  error,  an  addi¬ 
tional  option  is  provided  for  up  to  n  corrections  and  re- corrections, 
depending  upon  the  magnitude  of  the  quaternion  error,  €  ,  as  defined  in 

the  previous  section.  If  €  >  additional  corrections  and  re¬ 

corrections  will  be  made  until  the  number  of  corrections  (NCOR)  equals 
the  specified  maximum  number  of  corrections  (NCMAX). 

The  accumulative  integration  error  for  selected  variables  is 
provided  in  the  optional  program  output,  based  on  the  errors  computed  in 
the  last  correction  or  re- correction. 


where 

W 

(Xi\ 


SAVE 


L.  OPTIONAL  APPROXIMATE  EQUATIONS  OF  MOTION  FOR 
USE  WHEN  THE  NUTATION  IS  DAMPED 

Approximate  equations  of  motion,  with  q  =  r  =  0 ,  may  be 
selected  at  a  particular  time,  if  the  fixed-plane  axes  option  is  exercised. 

In  effect,  this  operation  eliminates  the  high  frequency  oscillatory  motion 
about  the  lateral  axes,  which  is  usually  associated  with  the  nutational  mode. 
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The  angular  velocities  q  and  r  are  determined  by  solution  of  the  differ¬ 
ential  equations  which  result  with  q  =  r  =  0  .  These  modified  equations 
are : 


The  minus  sign  is  selected  for  the  radical,  corresponding  to  the  slow  pre- 
cessional  mode.  The  positive  root  represents  a  flat  spin  mode  with  large 
angular  rates  and  is  of  little  interest. 

•  * 

The  optional  equations  of  motion  for  q  =  r  =  0  are  initiated 
at  a  time  when  the  nutational  oscillations  are  known  to  be  damped  to  an 
amplitude  of  a  few  degrees.  The  resulting  motion  predictions  are  valid  to 
the  extent  that  the  nutational  motion  can  be  neglected.  The  approximations 
are  invalid  if  the  high  frequency  motion  is  undampted. 

Because  the  angular  rates  will  usually  be  small  with 
q  =  r  =  0,  a  larger  integration  interval  can  be  utilized,  and  an  option  for 
this  new  time  interval  is  included  in  the  program  input. 

M.  AUXILIARY  FUNCTIONS 

For  purposes  of  program  output,  the  Euler  angle  rates  are 
specified,  although  these  quantities  do  not  enter  into  the  equations  of  motion. 
The  following  relationships  are  utilized  for  computing  the  Euler  angle  time 
derivatives'. 
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Description  of  Labeled  Common  Names 
Description 


Name 

(bTOf)  Arguments  of  SETMAX 

XMAT(3,3)  Rotational  Matrix 


Reference 
Section  -  Page 


^MPGIN^  Input  Related  Only  to  Main  Program 


XDOT(3) 

X.Y.Z 

II 

31 

P 

P 

II 

31 

Q 

q 

II 

31 

R  ' 

r 

n 

31 

ZALT 

Z 

H 

31 

PHI 

<P 

H 

31 

THETA 

0 

II 

31 

PSI 

r 

H 

3i 

TSTEP 

At 

n 

32 

INCPT1 

Print  cycle  before  q  =  »  =  0 

II 

32 

TMAX 

lmax 

II 

32 

ZSTOP 

^min 

H 

32 

TCH 

'"change  (q  -  T  -  0) 

II 

32 

TNEW 

At2 

H 

32 

NRK 

Integration  code 

II 

29 

NCMAX 

Number  of  corrections 

II 

29 

EPSMAX 

£ 

cmax 

H 

29 

ZETD1 

$,»  deg. 

II 

32 

ZETD2 

deg. 

II 

32 

INCPT2 

Print  cycle  after  q  -  *  =  0 

H 

32 
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description  of  Labeled  Common  Names  (Continued) 


Name  Description  Reference 

Section  -  Page 


<  allin)  Input  Related  to  Both  Main  Program  and 

Equations  of  Motion 

DIX 

*x 

II 

30 

DI 

1  or  Iy 

II 

30 

DIZ 

I* 

II 

30 

DIXY 

w 

11 

30 

BODFIX 

T/F  Flag  False  if  Fixed  Plane 

II 

30 

^EQMPG^  Main  Program  Constants  Needed  In 

Equations  of  Motion 

OR 

T/F  Flag  -  if  false,  use  simplified 

equations  for  q  and  r 

RAT 

W* y 

I 

9 

BODJX 

*xyAx 

I 

8 

BODJY 

*sy/xy 

I 

8 

bowz  i^nz 

I 

8 

BODEN 

1  -  *Xy^/ 

I 

8 

BODJ  <IX  +  Iy  -  Iz)  Ixy/dxly) 

i 

8 

-.BODPD 

^xy  ^y  {^y  ” 

/dxly) 

i 

8 

BODQD 

*xy  ~  *x  (*x  “  I a ) 

i 

8 

BODRD  (lx  -  Iy)/Iz 

i 

8 

ZET1 

S,  radians 

i 

21 

ZET2 

radians 

i 

21 

— -e8criPtlon  of  J-abel3d  Common  Names  (Continued V 

Name  n ^  Reference 

Ocscnphon  Section  -  Pa*. 


<  EQMIN>  Input  Related  Only  to  Equations  of  Motion 


IMAX 
JMAX 
KMAX 
LMAX 
TABA  (20) 
tabM(5) 
TABAR  (5) 
TABZ  (10) 

CX  (20,  5) 

CN  :(20,  5) 

CM  (20.  5) 
CMQ  (20,  5) 
CNR  (20,  5) 
CMPR  (20,  5) 
CNPQ  (20,  5) 
CNPA  (20,  5,5) 
CL  (20;  5,  5) 
CLP  (20,  5,5) 
CMPA  (20,  5,  5) 
CYO  (20,  5) 
CZO  (20,  5) 
CMO  (20,  5) 
CNO  (20,  5) 
CSF1  (20,  5) 
CN1  (20,  5) 
CSF3  (20~5) 


—  cv 

II 

29 

£5 

II 

29 

<5 

II 

29 

±  10 

II 

29 

(<*-) 

II 

35 

(m) 

II 

35 

(pd/2v) 

n 

35 

(-  z) 

n 

41 

cx 

n 

36 

CN 

n 

36 

CM 

n 

36 

Cmq 

ii 

3? 

cnr 

n 

37 

CmPr 

n 

38 

CnPq 

n 

38 

CNp  =  -CLp 

n 

39 

C* 

u 

39 

C  6 

*p 

n 

40 

p 

n 

40 

C7o 

n 

42 

o 

N 

n 

42 

Cm0 

n 

43 

C«o 

n 

43 

cSFj 

n 

44 

cNi 

n 

44 

cSF3 

ii 

44 

54 


Name 


Description  of  Labeled  Common  Names  (Continued) 

Reference 

Description  Section  -  Page 


{  EQMIN^  (continued) 


CSM1  (20,  5) 

cSMi 

H 

45 

CM1  (20,  5) 

CMl 

II 

45 

CSM3  (20,  5) 

cSM3 

II 

45 

CLPH1  (20,  5) 

CUl 

II 

46 

CLPH2  (20,  5) 

cUz 

H 

46 

WDX  (10) 

xw 

H 

41 

WDY  (10) 

Yw 

H 

41 

G 

g 

H 

30 

DMM 

m 

II 

30 

S 

S 

H 

31 

DEE 

d 

H 

31 

DY 

Ay 

II 

32 

ETA1 

Hi 

II 

32 

ETA2 

^2 

H 

32 

(  OUTIN )  Input  Related  Only  to  Output 

HEADER.  (11) 

IRUN 

IDATE 


H  30 


30 

in 


TT 


Name 


Description  of  Labeled  Common  Names  (Continued) 

Reference 

Description  Section  -  Page 


^TBL^  Arguments  for  TBL1,  TBL2,  TBL3 
found  from  GE TUN 


N1 

nl 

index 

such  that  . 

XT  (n 

-1)< 

XT  (n) 

N2  ; 

n2 

for  n 

=  nj  or  ng 

or  n3 

N3 

n3 

PI 

1 

ratio 

of  xT.feL 

°  XT  (n) 

-  X 

-  XT  (n 

P2 

2 

P3 

3 

for  n 

=  nj  or  n2 

or  n3 

(  OUTEM  V  Equations  Of  Motion  Generated  Data  for  Output 


CAPHI1 

*4-(S-s.) 

I 

16 

CAPHI2 

V  VCs-sJ 

I 

16 

PDV 

pd/2V 

I 

10 

ALPHA 

ot. 

I 

10 

CAPVA 

VA 

I 

20 

EM 

Mach  Number 

I 

10 

PEE 

p  or^-r  tan 8") 

I 

7,  25 

56 


Hicrachy  of  Programs  and  their  Labeled  Common 


Key 


Program  Names] 


(Arguments) 


^Labeled  Common \ 
'  Names  / 


•4rEPRT  (M£OF)| 

-»|readin1 

<MPGlN,  ALLIN,  EQMIN,  OUTIN^ 


•>|runge  (D,  F,  DT,  G)| 


SIXDEG 


<BTOF,MPGIN, 
ALLIN,  EQMPG> 


V — t 

\ 

’ 

;-*|eQMOT 

(Y,Z)h - t — — » 

\  ALLIN,  EQMPG, 
EQMIN,  OUTEM,  TBL> 


SETMAX(XL)| 

<BTOF> 

[—^NORMAL  (XL)  | 


GETLIN  (X.XTAB, 
NMAX.  N.  P) 


— H  TBL  1  (TAB)  I 
<TBL> 


— »l  TBL  2  (TAB)| 


<TBI> 

1 — >TtbL  3  (TAB)l 


<TBL> 


HWRITIT  (TIME,  Y,  E,  Z.  EPSIL,  LINE)  j 
<BTOF,  OUTIN,  OUTEM> 
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Flow  Chart  of  Main  Program 


if  Fixed- 
Plane  Axes 


59 


iVlain  Program  (continued) 


(100? 

J  Calculate  Equations  of  Motion 
1  F  (i,  NST) 


if  -  Z  >  Zg-^OP  °r  t  k  t; 


max 


if  continue 
and  no  print  I 


if  continue 
I  and  print 


Print  Output 


Print  Output 


<D 


for  new  run 


if  time  to  change  to 

steady?  state  calculations  for  q,  r 


if  not 


Set  Flag  for  q  r  False 
Set  At  to  new  &t 
-Set-print  index-to  2nd  index  I 


Bump  time,  print  index 
and  q  r  change  index 


Calculate  Runge-Kutta 
value  for  new  time 
D(i,  NST  +  1) 


if  NST  =  4  (i.  e. : 

Runge-Kutta 

finished) 


D(i,  4)  -*D(i,  Nstart)| 

Ei  =  0  _  '  __ . 


|  Calculate -:Equatiohs  of  Motion-; 
F(i,  Nstart)  to  clear  1st 
iteration 


for 

restart 


NST  =  NST  +  1 


if  NST  <  4 
or  Trapezoid| 


-M200 


for 

prediction 


NST  =  4 

and  Adams  or  Milne 


Calculate  Runge-Kutta  value 
for  new  time  from  initial 
time  to  estimate  E^  error 
D(i,5) 


for 

-correction 


60 


Flow  Chart  of  READIN 


Flow  Chart  of  RUNGE 


Flow  Chart  ofEQMOT 


mm 

Z  = 

u 

• 

u 

V 

• 

V 

W 

• 

w 

P 

• 

p 

q 

• 

q 

r 

• 

r 

X 

• 

X 

Y 

• 

Y 

Z 

• 

Z 

*0 

*0 

Xx 

*2 

x2 

*3 

x3  * 

«* 


c 


RE  TU  RN 
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EQMOT  Equations  in  More  Detail  for  Fx,  Fv>  Fz,  L,  M,  N 

denotes  body-fixed  axes  option  only 


- - 


CON  =  \  p  V2S 

Fx  =  CON*  Cx  (  oc  ,  M) 

Cm  ADD  =  .[+]  CNp  (  M  »  fy  2V 

r 


+  CSFl  (  &  ,  M)  sin  (Yl!  [f -S  +  Si3)  +  CSFs  (  5. ,  M) 

f ,  com  |  Ic'Js.'m)]  -  jcNia,  M)  +  1 

lin  j\ 

.  CON*{  -[cNta.  M)  ♦ 

-  c^p  ADD  •  cos  §\ 


+  Cn^ADD  *  sin  5  >- 


L  =  CON  d* 


+  |ci{  raTMVsi'inrif-S^S  +  ci,.  (S.M)  »in 


/  Cj  (&,  M.  ||)  +  pf‘  Cfj 


<&.  M,f|) 


I - 1 

-jAy  FZ  | 

U  - » 
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Flow  Chart  of  WRITIT 


1 


9 , 0  ,  Y  ,  X 

^of  errors,  Ap,  Aq,  Ar,  A^  ; 
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-r  low  i^nart  ot  SKTMAX 


SETMAX  (XL)  • 

_  * _ _ 

Given  XL,  —XL4  (  A,  — »A.) 

XMAT  (1,  1)  =  Xj  tXj  -Xj  Xj 

xmat  (2,  2)  =  x;  vXf  +x;  -x; 

XMAT  (3,  3)  =  A!  -  X,*  -XI  4-Xj 
XMAT(I,  2)  =  2  (  X,Xt  -  X.X,) 
XMAT  (2,  1)  =  2(  X,  X»  +  X,X,j 
XMAT.(i,  3)  ,  2  (  A,  A,  +  X.X.) 

XMAT  (3.  i)  =  2  (  X,Xj  -  A.  A.) 
XMAT  (2,  3)  =  2(  X.X,  -  X.X,) 
XMAT  (3.  2)  e  2  (  X,X,  3.  X.X,) 


^  return 
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F  low 


Ch  a,r  t 


of  GETLIN 


Note:  N  &  P  are  used  such  that 

Y  -  (1  -  P)«Yn  +  P  •  Yn  _  x  in  TBL.  1,  2  and  3 

*  i.  e; ,  to  handle  negative  altitude  and  negative  p  where  only  +  arguments 
are  in  Table. 

**  i.  e, ,  saves  interpolation.  *  - 
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Flow  Chart  of  TBL  1,  TBL  2,  TBL  3 


*  TBL  2  =  (1-P2)  t(l-Pl).  TAB  (N1.N2)  +  PI-  TAB  (Nl-1,  N2)J 

+  P2[(lrPl)-TAB  (N1.N2-1)  +  Pl'TAB  (Nl-1,  N2-1)J 

**  TBL  3  =  (l-P3)  Kl-P2){(l-Pl)* **TAB  (N1.N2.N3)  +  P1-TAB(N1- 1,  N2,  N3)} 

+  P2{(l-Pl)*TAB  (Nl,  N2- 1,  N3)  +  PI 'TAB  (N1-1.N2-1.N3)}} 

.  +  P3  [(1-PZ)^  (l-Pl)'TAB  (N1.N2.N3-1)  +  PI-TAB  (Nl- 1,  N2,  N3- 1)} 
+  P2{(1-P1)»TAB  (Nl,  N2r  1,  N3-  1)  +  Pl'TAB  (N2- 1,  N2- 1,  N3- 1)}] 
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Compact  Input  Sheet 


No.  of 

N  Columns 


- 1 

3511 

KRD(I)  I-  1,  35 

1 

35 

0 

13, 

2.X11A6, 

I  Run,  HEADER,  KA6,  BOD  FIX-,  I 

Lead  Card  1 

78 

1X11,  LI, 

14 

IDATE 

1 

1 

612, 

F12.  5 

IMAX,  JMAX,  KMAX,  LMAX,  ! 

Subscript  Card 

24 

NRK,  NCMAX,  EPSMAX 

2 

6F12..5 

TABA(I) 

12- 

72 

3 

TABM(J) 

4 

TABAR(K) 

5 

CX  (I,  J)  J  then  I 

6 

CN  (I,  J)  J  then  1 

7 

CM  (I,  J)  J  then  I 

8 

CMQ  (I,  J)  J  then  I 

9 

CNR  (I,  J)  J  then  I 

10 

CMPR  (I,  J)  J  then  I 

11 

CMPQ  (I,J)J  then  I" 

12 

CNPA  (I,  J,  K.)  K  then  J 

then  I 

* 

13 

CL  (I,  J,  K)  K  then  J 

then  1 

14 

CLP  (I;  J,  K)  K  then  J 

then  I 

15 

CMPA  (I,  J,  K)  K  then  J 

then  I 

16 

G,  DMM,  DIX,  :DI,  DIZ,  DIXY 

j  Lead  Card  2 

72 

17 

S,  DEE,  THETA,  PSI,  Z,  PHI 

J  Lead  Card  3 

72 

18 

XDOT,  YDOT,  ZDOT,  P,  Q,  R 

j  Lead  Card  4 

72 

19 

F12.5,  4X18, 

TSTEP,  INCPT1,  TMAX,  ZSTOP, 

j  Lead  Card -5 

76 

4F12.5,  14 

TCH,  TNEW„,  INCPT2 

20 

6F  12;5 

TABZ(L). 

12- 

72 

WDX(L) 

12- 

72 

WDY(L) 

_ =._ 

12- 

72 

21 

6F12.5  J 

i' 

DY.ZETDl,  ZETD2 

|  Lead  Card  6 

60 

ETA1.ETA2 

• 

22 

“ 

' 

CYO  (I,  J)  J  then  I 

12- 

72 

23. 

CZO  .(I,J)  J  then  I 

24 

Body 

CMO  (I,  J)  J  then  I 

25 

Axis 

CNO  (I,  J)  J  then  I 

26 

Only 

CSFI  Cl,  J )  J  then  I 

27 

CN1  U.  J)  J  then  I 

28 

CSF3  (I.  J)  J  then  I 

29 

CSMl  (I.  J)  J  then  I 

3C 

CM1  <i»  J)  J  then  I 

31 

CSM3  (I,  J)  J  then  I 

32 

CLPM1  (l.J)  J  then  I 

33 

' 

f  \ 

t 

CLPM2  (I.  J)  J  then  I 

—  '  -  * -  - -  * 

N  is  a  code  for  card  type  (helpful  to  punch  in  columns  79-80) 

N  >  6  all  optional  0.1  added  funs  (KRD(N)  =  0  no,  KRD(N)  =  1  read) 

N  =  20  also  optional  on  any  fun  (LMAX  =•  0  no,  LMAX  >  6  read) 

N  =  0  required  for  each  run. 

N  =  - 1  required  for  each  additional  run. 

Cards  0  —*•19  required  for  FIXED  PLANE  1st  run;  20  if  LMAX  >  0. 
Cards  0-*19,  21-33  required  for  BODY  FIXED  1st  run;  20  if  LMAX  >  0. 
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Z<11>  =  .5*<  Y(iO)*PEE-Y(13)*t13)  +  t(i?)*Y(6)  )  EQMOT 
Z(12)  =  .5*(  Y(l3)*PFF  +  Y(in)*Y(5>-Y(H)#Y(6))  EQMOT 
Z(l3)  =  .5*(-.Y(i2)*PEE+Y(ll)*Y(5)  +  y<iO)*Y<6)  )  EQMOT 
RETURN  EQMOT 
END  EQMOT 
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SECTION  VI 


COMMENTS  AND  SPECIAL  INSTRUCTIONS 


A.  MAGNUS  ROTOR  TRAJECTORY  AND  MOTION  SIMULATION 

In  most  instances,  magnus  rotor  trajectories  will  be  computed 
with  fixed-plane  axes.  Where  complete  trajectory  data  are  to  be  obtained 
(launch  to  impact)  for  ilight  times  of  several  seconds  or  more,  it  will  often 
be  convenient  to  use  the  q  =  r  =  0  option  after  about  the  first  second  of 
flight*  (or  after  a  time  where  the  nutation  is  damped  to  a  few  degrees  am¬ 
plitude).  For  the  investigation  of  auto  rotation  initiation  and  the  effects  of 
dynamic  unbalance,  the  body-fixed  axes  may  be  used. 

For  initiating  magnu?  rotor  trajectories,  it  is  convenient  to 
use  the  initial  horizontal  velocity  component  along  the  (  -Y )  axis,  such  that 
the  Euler  angles,  9  and  Y  ,  will  approach  zero  if  the  rotor  is  in  gliding 
flight  attitude,  3L  -*>  'fr/2  .  In  this  manner,  positive  spin  rates,  correspond¬ 
ing  to  positive  values  of  the  spin  torque  coefficient  Cj^(  $.  ,  M,  ££ ),  will 
result  in  an  upward  magnus  lift  force. 

Special  care  must  be  used  in  selecting  the  initial  conditions. 
Letting  the  initial  flight  path  angle,  Y  ,  be  specified  by 

X  =  0 

-  Y  =  V  cos  X  1) 

z  =  v  sin  r 


then  the  initial  Euler  angles,  9  and  Y  *or  the  fixed-plane  axes,  are 
related  to  the  initial  sideslip  angle,  ,  and  the  orientation  of  the  angle  of 
attack  plane,  ,  by. the  following  relations  (see  also  Figure  6). 

sin  (  -  0  )  =  cos^  sin  <p  cos  jl*  +  sin^  sin)4  2) 

,  ,  .  -cos  /3  sin  <f  sin)"  +  sin,a  cos  Y 


tan  (  -  y  )  - 


cos 0  cos 


where 

-j  <  Y  c 

£ 

Z 

for 

-f< 

<5 

and 

^  <  Y  < 

£ 

z 

for 

£  < 

2 

*  Note:  if  the  q  =  r  =  0  option  is  not  used  with  fixed  plane  axes,  then 
TCH  >  TSTOP  for  input  lead  card  5. 
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Figure  6.  Description  of  Initial  Conditions  as  a  Function  of  X ,  /3  ,  <P  and  A 


If  an  initial  coning  motion  is  desired,  such  that  the  total 
angular  momentum  does  not  coincide  with  x  axis,  then  the  coning  angle 
can  be  specified  by,  X  «  and  the  equivalent  cross-angular  velocity  deter¬ 


mined  by 


4) 


The  angular  velocity  components  with  respect  to  the  lateral 
fixed-plane  axes  are  given  in  terms  of  Slan  by  the  relationships 

vvj* 


q  =  jfl, A  cos  -  A.)  [cos$  sin  (-  r)  +  sin*  sin  /  cos  (-V)] 

eq- » 


-sin  ()3  -  A  )  cos  /  cos  (  -  y ) 


eq.j  cos  (/?-*)  [-cos*  cos  ( 

-y)  sin  (-  0  ) 

+  sin  (j?  sin  /  sin  (~y  )  sin  {-  Q  )  +  sin  *  cos  /  cos  (-  0  )] 

+  sin  (/S-A)[  -cos  /  sin  (-  y)  sin  (-  0  )  +  sin  /  cos  (-  0  )] 


When  random  values  of  <2  and  $>  are  to  be  used  for  Monte 
Carlo  simulation  of  impact  patterns,  relations  2)  and  3)  will  be  used  exten¬ 
sively.  In  simulation  of  impact  patterns,  symmetry  considerations  can 
also  be  employed  such  that  only  the  left  or  right  hand  side  of  the  pattern 
need  be  determined.  The  symmetry  considerations  are  applied  to  the 
initial  conditions  by  restricting  the  angle  of  attack  to  O  <oc  Z  j-  . 

B.  BODY-FIXED  AXES 

Ballistic  trajectories  of  slowly- spinning  rockets,  missiles, 
bombs,  and  projectiles  will  usually  be  computed  using  the  body-fixed  axes 
option.  For  ballistic-type  trajectories  it  will  be  convenient  to  select  the 
horizontal  component  of  the  initial  velocity  in  the  direction  of  the  positive 
X  axis,  in  order  that  first  quadrant  values  can  be  used  for  the  Euler  angles, 
0  and  r. 

Supplementary  calculations  will  be  required  to  determine  the 
initial  Euler  angles  if  only  the  velocity  vector  and  the  angle  of  attack  param¬ 
eters,  <2  and  $>  are  known  or  specified.  The  procedure  involves  finding 
the  values  of  0  and  *y  which  satisfy  the  matrix  equation 
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COSY 


-so 


r  i 

u 


w 


co-cy  cosy  - so 

CY’Sy -so-ay-sY  sp  so-sr+  cqcy  sq>c$ 
cy-cq-so  +  sO'SY  CQ-so-sY-sfp'Cf  c<p>co 


subject  to  the  constraints 


oC 

$ 


=  tan 


-l  J/ra+  url 


u 


-  tan 


1  ar 


fc  = 


jr 

CONSTANT 


X 

Y 

Z 


For  special  cases,  closed  form  solutions  exist,  i.  e. ,  for  Y  = 


Z  =  <Pe  =  0 


tan  Q  = 


tarn,  oi 


COS  *Y  - 


co<-  6 
Ctt<-  8L 


More  general  velocity  vector  orientations  with  4  0  require  complicated 
numerical  solutions.  Obviously,  these  problems  do  not  exist  if  the  Euler 
angles  can  be  selected  apriori.  However,  for  multiple  simulations,  involv¬ 
ing  a  fixed  velocity  vector,  the  above  procedure  will  normally  be  required. 
A  nearly  general  solution  for  0  >Y  >  subject  only  to  Y  =  0,  i.  e. , 
tan  /  =  Z/X,  is  given  by  the  relations 

sin  6  [-  sin  t  J  +  cos  6  [  cos  t  cos  Yl  •  es  cos  ^ 


sin  0 


f  j  1  +  B  f  tan  y  I  _  sin  4>.  -f  cos  <P. 
L  tan  Iff  j  [  sin  Y  J  sin  tan  S  “ 


tan  5 

COS 


100 


c. 


INTEGRATION  INTERVAL 


For  magnus  rotor  motion  simulations,  it  will  generally  be 
advisable  to  select  an  integration  time  interval  (TSTEP)  such  that  about  100 
integrations  are  performed  for  each  cyclic  period.  For  fixed-plane  axes, 
the  highest  frequency  will  be  the  nutation  frequency,  while  for  body -fixed 
axes  the  highest  frequency  will  correspond  to  the  spin  rate.  Since  the  nuta¬ 
tion  frequency  is  approximately  plx/l,  a  large  saving  in  computing  time 
can  be  achieved  using  the  fixed-plane  axes  if  Ix/I  «  1  . 

When  it  is  anticipated  that  only  a  short  segment  of  a  motion 
history  will  involve  high  frequencies,  localized  integration  improvement  can 
be  achieved  by  using  the  re-correction  option.  A  quaternion  error  param¬ 
eter  EPSMAX  of  about  1.0  x  10"  would  then  be  used  in  conjunction  with 
a  value  of  NCMAX  of  3  or  4. 
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control  inputs,  aerodynamic  roll  angle  effects,  flow  asymmetry  with  respect  to  the 
angle  of  attack  plane  at  zero  spin,  and  lateral  c.  g.  offset.  ~-5 

-  The  computer  program  is  written  in  General  Fortran  IV  langu  age  compatible 
with  CDC  6400,  IBM  360,  and  GE  635  data  processing  machines.  Included  in  the 
report  are  the  program  input  and  output  formats,  flow  charts  of  the  main  program 
subroutines,  and  a  complete  program  listing.  / 
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